using System;
using DambachMulti;
using Real = System.Double;

namespace LinearAlgebra.Matrices
{
	public partial class Matrix
	{
		/// <summary>
		/// Set every element in a matrix to a value given by (i,j)=>{ ... }
		/// </summary>
		/// <param name="f">The function</param>
		public Matrix Set( MatrixElementFunction f )
		{
			for (int i = 0; i < this.Length; i++)			
				for (int j = 0; j < this[0].Length; j++)				
					this[i][j] = f( i, j );			
			return this;
		}

		/// <summary>
		/// Let each element (i,j) be defined by a function f(x) given by the user
		/// </summary>				
		public Matrix Set( RealFunction f )
		{
			for (int i = 0; i < this.Length; i++)			
				for (int j = 0; j < this[0].Length; j++)				
					this[i][j] = f( this[i][j] );			
			return this;
		}

		/// <summary>
		/// Do something for each coordinate of a matrix
		/// </summary>		
		public Matrix Do( MatrixElementProcedure proc )
		{
			for (int i = 0; i < this.Length; i++)			
				for (int j = 0; j < this[0].Length; j++)				
					proc( i, j );
			return this;
		}
		
		/// <summary>
		/// Do something for each value of a matrix
		/// </summary>		
		public Matrix Do( RealProcedure proc )
		{
			for (int i = 0; i < this.Length; i++)			
				for (int j = 0; j < this[0].Length; j++)				
					proc( this[i][j] );
			return this;
		}		

		/// <summary>
		/// Add a matrix of the same row height to the end of this matrix
		/// </summary>
		/// <param name="B">The matrix to append to this one</param>
		/// <returns>A reference to the modified matrix</returns>
		public Matrix Augment( Matrix B )
		{
			if( this.RowCount != B.RowCount )
				return this;
			for( int i = 0; i < this.RowCount; i++ ) 	//for every row
				for ( int j = 0; j < B.ColumnCount; j++ )
					this[i].Add( B[i][j] );
			return this;
		}

		/// <summary>
		/// Returns a new matrix that has had the augmentation removed
		/// </summary>
		/// <param name="n"></param>
		/// <param name="m"></param>
		/// <returns></returns>
		public Matrix DeAugment( int m )
		{
			Matrix result = GetSubmatrix( 0, this[0].Length-m, m, m );
			for( int i = 0; i < this.RowCount; i++ )
				this[i].RemoveRange( this[i].Length-m, m );			
			return result;
		}

		/// <summary>
		/// Creates a new copy of the matrix
		/// </summary>
		/// <returns></returns>
		public Matrix Copy()
		{
			return this.GetSubmatrix( 0, 0, RowCount, ColumnCount );
		}
		
		/// <summary>
		/// Think substring, only in 2-dimensions instead of one
		/// </summary>
		/// <param name="startingRow">The starting row (inclusive)</param>
		/// <param name="startingCol">The starting column (inclusive) </param>
		/// <param name="rowCount">The number of rows from the starting row to include.</param>
		/// <param name="colCount">The number of columns from the starting column to include</param>		
		public Matrix GetSubmatrix( int startingRow, int startingCol, int rowCount, int colCount )
		{
			Matrix result = new Matrix( rowCount, 0 );
			this.Do( delegate( int i, int j ) { 
				if( i >= startingRow && 
					j >= startingCol && 
					i < (startingRow + rowCount) && 
					j < (startingCol + colCount) ) {
					result[i].Add( this[i][j] );
				}
			} );
			return result;
		}
		
		/// <summary>
		/// If this matrix is invertible, this function returns it's inverse, otherwise it throws an exception
		/// </summary>		
		public Matrix GetInverse()
		{
			return GetInverse( Real.Epsilon );
		}

		/// <summary>
		/// If this matrix is invertible, this function returns it's inverse, otherwise it throws an exception
		/// </summary>		
		public Matrix GetInverse( Real epsilon )
		{			
			Matrix copy = this.Copy();
			int n = this.RowCount; //only square matrices are invertible

			copy.Augment( Identity( n, n ) );
			copy.Eliminate( epsilon );
			copy.Normalize(); //set all pivots to 1 and maintain row equivalence
			copy.SolveByBackSubstitution();
			return copy.DeAugment( n );
		}
		
		/// <summary>
		/// Takes this matrix to Row Echelon Form
		/// </summary>				
		public Matrix Eliminate() { return Eliminate( Real.Epsilon ); }

		/// <summary>
		/// Performs standard Gaussian Elimination to bring the matrix to Row Reduced-Upper Triangular form
		/// </summary>
		/// <param name="zeroEpsilon">If a pivot value is smaller than the zero epsilon, a divide by zero exception is thrown</param>
		/// <returns>A Row Reduced Matrix</returns>
		public Matrix Eliminate( Real zeroEpsilon )
		{
			for( int k = 0; k < this.Length; k++ ) //for every row
			{
				// Check for small element on the diagonal
				if ( Math.Abs( this[k][k]) < zeroEpsilon ) //C# has this built in exception type 
					throw new DivideByZeroException( "Value less than epsilon(" + zeroEpsilon + ") on the diagonal, possible divide by zero condition! Cannot continue." );
				
				Multi.For( k + 1, this.Length, (i)=> {  //Get zeros for column k, 												
					//Real xmult; //WARNING: efficiency optimization
					//recurse over every row under the kth row
					Real xmult = -this[i][k] / this[k][k];	//what is the mult?
					this.AddRowMultiple( xmult, k, i );	//reduce mat this, inverse mat B
				} );
			}
			return this;
		}

		/// <summary>
		/// Produces a matrix that is in reduced row echelon form (row canonical form)
		/// </summary>		
		public Matrix RowReduce() {
			return RowReduce( Double.Epsilon );
		}

		/// <summary>
		/// Produces a matrix that is in reduced row echelon form (row canonical form)
		/// </summary>		
		public Matrix RowReduce( Real epsilon ) { return Eliminate( epsilon ).Normalize().SolveByBackSubstitution(); }
		
		/// <summary>
		/// Divides every element on a row by the pivot value for that row
		/// </summary>
		/// <returns>A reference to the normalized matrix</returns>
		public Matrix Normalize()
		{
			Real piv; //this is to reduce the garbage collection = faster code			
			for( int i = 0; i < this.RowCount; i++ ) { //this loop performs normalization of the matrix...
				piv = this[i][i];
				for( int j = 0; j < this.ColumnCount; j++ ) {
					this[i][j] /= piv;
				}
			}
			return this;
		}
		
		/// <summary>
		/// Takes a Matrix that is in Row Reduced Form and Solves it
		/// </summary>
		/// <returns></returns>
		public Matrix SolveByBackSubstitution()
		{			
			for( int j = this.RowCount-1; j >= 0; j-- )	//go backwards and solve for other x values			
				for( int i = j-1; i >= 0; i-- )		//I don't think we can parallelize this part as is				
					AddRowMultiple( -this[i][j], j, i );
			return this;
		}
		
		/// <summary>
		/// Multiplies a Row Vector by a Scalar, and adds the result toRow the Target Row
		/// </summary>		
		public Matrix AddRowMultiple ( Real multiple, int fromRow, int toRow )
		{
			this[toRow] += (this[fromRow] * multiple);
			return this;			
		}
				
		public static Matrix operator*( Matrix A, Matrix B ) { return ParallelMultiply( A, B ); }

		/// <summary>
		/// Matrix Multiplication (Look it up dummy!)
		/// </summary>		
		private static Matrix ParallelMultiply( Matrix A, Matrix B )
		{
			if( A.ColumnCount != B.RowCount )
				throw new Exception("Matrices are of mismatched dimensions. Check that MxN * NxP == MxP");
			Matrix c = new Matrix( A.RowCount, B.ColumnCount );
			
			Multi.For( 0,  c.RowCount, 1, (i)=> {
				for( int j = 0; j < c.ColumnCount; j++ )
					for( int k = 0; k < A.ColumnCount; k++ ) 
						c[i][j] = c[i][j] + A[i][k] * B[k][j]; } );
			return c;
		}

		/// <summary>
		/// Matrix Multiplication (Look it up dummy!)
		/// </summary>		
		private static Matrix ReferenceMultiply( Matrix A, Matrix B )
		{
			if( A.ColumnCount != B.RowCount )
				throw new Exception("Matrices are of mismatched dimensions. Check that MxN * NxP == MxP");
			Matrix c = new Matrix( A.RowCount, B.ColumnCount );

			for( int i = 0; i < c.RowCount; i++ )
				for( int j = 0; j < c.ColumnCount; j++ )
					for( int k = 0; k < A.ColumnCount; k++ ) { 
						c[i][j] = c[i][j] + A[i][k] * B[k][j]; }
			return c;									
		}		

		/// <summary>
		/// Returns a matrix where each [i,j] = A[i,j] + B[i,j]
		/// </summary>
		public static Matrix operator+( Matrix A, Matrix B )
		{
			if( A.M != B.M && A.N != B.N )
				throw new Exception( "Matrices are not the same size, operation is undefined" );
			Matrix r = new Matrix( A.M, A.N );	
			r.Set( (i,j)=> A[i][j] + B[i][j] );
			return r;
		}
		
		/// <summary>
		/// Transposes the Matrix Elements in place
		/// </summary>		
		public Matrix GetTranspose() {			
			Matrix m = Create( this.N, this.M, (i, j)=> this[j][i] );
			return m;
		}				

		public Real GetDeterminant()
		{
			if( M != N )
				return 0;
			Matrix c = Copy().Eliminate();
			Console.WriteLine( c.ToString() );
			Real det = c[0][0];
			
			for( int i = 1; i < c.RowCount; i++ )			
				det *= c[i][i];
			return det;
		}
	}
}